Climate Water Loss Experiment - Treatment Analysis

Savannah Weaver

2022

Packages

if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse") # workflow and plots
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # for p-values
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr") # for multi-ggplot figs
if (!require("UsingR")) install.packages("UsingR")
library("UsingR") # simple.eda model assumption checker
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # lmer model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans') # marginal means
if (!require("AICcmodavg")) install.packages("AICcmodavg")
library("AICcmodavg") # model selection
if (!require("MuMIn")) install.packages("MuMIn")
library("MuMIn") # model selection
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # color
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown format

Background and Goals

This data was collected June - August by Master’s student Savannah Weaver, advisor Dr. Emily Taylor, and research assistants Tess McIntyre and Taylor Van Rossum. Adult male Sceloporus occidentalis were caught across the Cal Poly campus then acclimated to 4 different climate treatments. This R file analyzes the effect of experimental climate treatments on lizard body condition, osmotic balance (plasma osmolality and hematocrit), and cutaneous evaporative water loss (CEWL).

Data

Load

Read-in data that was compiled, formatted, and checked for completeness in ‘wrangling_general’. See that file for information related to the variables.

dat <- read_rds("./data/analysis_data_experiment.RDS")
summary(dat)
##  measurement_date        type     individual_ID     mass_g     
##  Min.   :2021-06-16   exp  :804   201    :  7   Min.   : 7.00  
##  1st Qu.:2021-07-01   rehab:132   202    :  7   1st Qu.: 9.50  
##  Median :2021-07-25               203    :  7   Median :10.60  
##  Mean   :2021-07-22               204    :  7   Mean   :10.64  
##  3rd Qu.:2021-08-14               205    :  7   3rd Qu.:11.70  
##  Max.   :2021-09-01               206    :  7   Max.   :17.40  
##                                   (Other):894                  
##  hematocrit_percent trial_number temp_tmt   humidity_tmt     SVL_mm     
##  Min.   :13.00      1:175        Hot :467   Humid:468    Min.   :60.00  
##  1st Qu.:26.00      2:203        Cool:469   Dry  :468    1st Qu.:66.00  
##  Median :32.00      3:231                                Median :67.00  
##  Mean   :31.99      4:189                                Mean   :67.74  
##  3rd Qu.:37.00      5:138                                3rd Qu.:70.00  
##  Max.   :52.00                                           Max.   :77.00  
##  NA's   :408                                                            
##          tmt          day_n        day_factor osmolality_mmol_kg_mean
##  Cool Humid:238   Min.   : 0.000   0 :134     Min.   :295.3          
##  Hot Humid :230   1st Qu.: 4.000   4 :134     1st Qu.:336.1          
##  Cool Dry  :231   Median : 6.000   5 :134     Median :351.3          
##  Hot Dry   :237   Mean   : 5.705   6 :134     Mean   :354.3          
##                   3rd Qu.: 8.000   7 :134     3rd Qu.:370.0          
##                   Max.   :10.000   8 :134     Max.   :471.5          
##                                    10:132     NA's   :414            
##  CEWL_g_m2h_mean   msmt_temp_C    msmt_RH_percent cloacal_temp_C 
##  Min.   : 7.152   Min.   :24.80   Min.   :25.52   Min.   :23.00  
##  1st Qu.:19.755   1st Qu.:26.20   1st Qu.:46.11   1st Qu.:25.00  
##  Median :24.152   Median :26.74   Median :47.88   Median :26.00  
##  Mean   :24.767   Mean   :26.72   Mean   :46.74   Mean   :25.92  
##  3rd Qu.:28.505   3rd Qu.:27.11   3rd Qu.:50.50   3rd Qu.:27.00  
##  Max.   :56.066   Max.   :29.20   Max.   :56.16   Max.   :30.00  
##  NA's   :669      NA's   :668     NA's   :668     NA's   :668    
##   msmt_temp_K      e_s_kPa_m       e_a_kPa_m       msmt_VPD_kPa  
##  Min.   :297.9   Min.   :3.219   Min.   :0.9894   Min.   :1.486  
##  1st Qu.:299.4   1st Qu.:3.504   1st Qu.:1.6464   1st Qu.:1.767  
##  Median :299.9   Median :3.620   Median :1.7411   Median :1.853  
##  Mean   :299.9   Mean   :3.620   Mean   :1.6833   Mean   :1.937  
##  3rd Qu.:300.3   3rd Qu.:3.701   3rd Qu.:1.7992   3rd Qu.:2.012  
##  Max.   :302.4   Max.   :4.194   Max.   :1.9326   Max.   :3.021  
##  NA's   :668     NA's   :668     NA's   :668      NA's   :668    
##       SMI         temp_mean_tmttrial temp_SD_tmttrial humidity_mean_tmttrial
##  Min.   : 6.747   Min.   :23.30      Min.   :0.5966   Min.   :13.75         
##  1st Qu.: 9.714   1st Qu.:24.05      1st Qu.:0.7828   1st Qu.:29.21         
##  Median :10.594   Median :24.88      Median :1.0461   Median :45.24         
##  Mean   :10.599   Mean   :29.60      Mean   :1.1513   Mean   :52.94         
##  3rd Qu.:11.390   3rd Qu.:35.05      3rd Qu.:1.5191   3rd Qu.:82.84         
##  Max.   :15.063   Max.   :36.00      Max.   :1.8447   Max.   :93.15         
##                                                                             
##  humidity_SD_tmttrial e_a_mean_tmttrial e_a_SD_tmttrial  VPD_mean_tmttrial
##  Min.   : 4.370       Min.   :0.4258    Min.   :0.1891   Min.   :0.2053   
##  1st Qu.: 6.234       1st Qu.:1.3896    1st Qu.:0.2651   1st Qu.:0.8018   
##  Median : 7.382       Median :2.2917    Median :0.3166   Median :2.0450   
##  Mean   : 8.765       Mean   :2.3340    Mean   :0.3936   Mean   :2.0043   
##  3rd Qu.:12.297       3rd Qu.:2.6749    3rd Qu.:0.5019   3rd Qu.:3.1651   
##  Max.   :19.846       Max.   :4.8058    Max.   :0.7691   Max.   :4.0685   
##                                                                           
##  VPD_SD_tmttrial      temp_C      humidity_percent    e_a_kPa    
##  Min.   :0.1690   Min.   :23.80   Min.   :17.90    Min.   :0.50  
##  1st Qu.:0.2558   1st Qu.:23.80   1st Qu.:34.10    1st Qu.:2.00  
##  Median :0.3483   Median :24.40   Median :56.20    Median :2.15  
##  Mean   :0.4038   Mean   :29.64   Mean   :52.84    Mean   :2.32  
##  3rd Qu.:0.4735   3rd Qu.:35.50   3rd Qu.:78.30    3rd Qu.:2.30  
##  Max.   :0.8863   Max.   :35.50   Max.   :80.90    Max.   :4.50  
##                                                                  
##     VPD_kPa     
##  Min.   :0.600  
##  1st Qu.:0.600  
##  Median :1.800  
##  Mean   :2.002  
##  3rd Qu.:3.800  
##  Max.   :3.800  
## 

Split

Make sub-dataframes without recovery data / with only recovery-related data:

note: recovery analysis was not included in publication

dat_no_rehab <- dat %>%
  dplyr::filter(day_n %in% c(seq(0,8))) 

recovery_values <- dat %>%
  dplyr::filter(day_n == 10) %>%
  dplyr::select(individual_ID,
                end_hct = hematocrit_percent,
                end_osml = osmolality_mmol_kg_mean,
                end_SMI = SMI)

recovery_v_post_exp <- dat %>%
  dplyr::filter(day_n == 8) %>%
  left_join(recovery_values, by = 'individual_ID') %>%
  mutate(delta_osml_10_8 = end_osml - osmolality_mmol_kg_mean,
         delta_hct_10_8 = end_hct - hematocrit_percent,
         delta_SMI_10_8 = end_SMI - SMI)

recovery_v_pre_exp <- dat %>%
  dplyr::filter(day_n == 0) %>%
  left_join(recovery_values, by = 'individual_ID') %>%
  mutate(delta_osml_10_0 = end_osml - osmolality_mmol_kg_mean,
         delta_hct_10_0 = end_hct - hematocrit_percent,
         delta_SMI_10_0 = end_SMI - SMI)

Check

Dates:

unique(dat$measurement_date)
##  [1] "2021-06-16" "2021-06-20" "2021-06-21" "2021-06-22" "2021-06-23"
##  [6] "2021-06-24" "2021-06-26" "2021-06-30" "2021-07-01" "2021-07-02"
## [11] "2021-07-03" "2021-07-04" "2021-07-06" "2021-07-20" "2021-07-24"
## [16] "2021-07-25" "2021-07-26" "2021-07-27" "2021-07-28" "2021-07-30"
## [21] "2021-08-08" "2021-08-12" "2021-08-13" "2021-08-14" "2021-08-15"
## [26] "2021-08-16" "2021-08-18" "2021-08-22" "2021-08-26" "2021-08-27"
## [31] "2021-08-28" "2021-08-29" "2021-08-30" "2021-09-01"

Number of measurements for each lizard:

dat_no_rehab %>%
  group_by(individual_ID) %>%
  summarise(n = n()) %>%
  arrange(n)
## # A tibble: 134 × 2
##    individual_ID     n
##    <fct>         <int>
##  1 201               6
##  2 202               6
##  3 203               6
##  4 204               6
##  5 205               6
##  6 206               6
##  7 207               6
##  8 208               6
##  9 209               6
## 10 210               6
## # … with 124 more rows

Every lizard has 6 experimental measurements: pre-tmt, mid-tmt, post-tmt, and mass checks on each of the 3 days between mid and post-tmt.

Did any of the treatment groups inherently start out with large differences in response variables?

dat %>%
  dplyr::filter(day_n == 0) %>%
  group_by(tmt) %>%
  summarise(mean(mass_g),
            sd(mass_g),
            mean(SMI),
            mean(hematocrit_percent),
            mean(osmolality_mmol_kg_mean),
            mean(CEWL_g_m2h_mean))
## # A tibble: 4 × 7
##   tmt        `mean(mass_g)` `sd(mass_g)` `mean(SMI)` mean(hema…¹ mean(…² mean(…³
##   <fct>               <dbl>        <dbl>       <dbl>       <dbl>   <dbl>   <dbl>
## 1 Cool Humid           11.6         1.35        11.7        39.6    351.    20.9
## 2 Hot Humid            11.6         1.75        11.5        37.9    347.    21.4
## 3 Cool Dry             11.8         1.61        11.8        39.3    346.    20.0
## 4 Hot Dry              12.0         1.68        11.8        38.9    347.    20.9
## # … with abbreviated variable names ¹​`mean(hematocrit_percent)`,
## #   ²​`mean(osmolality_mmol_kg_mean)`, ³​`mean(CEWL_g_m2h_mean)`

There are slight differences, but overall the starting values across groups are more or less the same.

Temp & RH during (all, before and after exp) CEWL measurements:

summary(dat_no_rehab)
##  measurement_date        type     individual_ID     mass_g     
##  Min.   :2021-06-16   exp  :804   201    :  6   Min.   : 7.00  
##  1st Qu.:2021-06-30   rehab:  0   202    :  6   1st Qu.: 9.50  
##  Median :2021-07-25               203    :  6   Median :10.60  
##  Mean   :2021-07-22               204    :  6   Mean   :10.65  
##  3rd Qu.:2021-08-13               205    :  6   3rd Qu.:11.60  
##  Max.   :2021-08-30               206    :  6   Max.   :17.40  
##                                   (Other):768                  
##  hematocrit_percent trial_number temp_tmt   humidity_tmt     SVL_mm     
##  Min.   :13.00      1:150        Hot :402   Humid:402    Min.   :60.00  
##  1st Qu.:28.25      2:174        Cool:402   Dry  :402    1st Qu.:66.00  
##  Median :33.00      3:198                                Median :67.00  
##  Mean   :33.75      4:162                                Mean   :67.73  
##  3rd Qu.:39.00      5:120                                3rd Qu.:70.00  
##  Max.   :52.00                                           Max.   :77.00  
##  NA's   :406                                                            
##          tmt          day_n     day_factor osmolality_mmol_kg_mean
##  Cool Humid:204   Min.   :0.0   0 :134     Min.   :295.3          
##  Hot Humid :198   1st Qu.:4.0   4 :134     1st Qu.:334.7          
##  Cool Dry  :198   Median :5.5   5 :134     Median :348.3          
##  Hot Dry   :204   Mean   :5.0   6 :134     Mean   :352.3          
##                   3rd Qu.:7.0   7 :134     3rd Qu.:367.4          
##                   Max.   :8.0   8 :134     Max.   :445.5          
##                                 10:  0     NA's   :413            
##  CEWL_g_m2h_mean   msmt_temp_C    msmt_RH_percent cloacal_temp_C 
##  Min.   : 7.152   Min.   :24.80   Min.   :25.52   Min.   :23.00  
##  1st Qu.:19.755   1st Qu.:26.20   1st Qu.:46.11   1st Qu.:25.00  
##  Median :24.152   Median :26.74   Median :47.88   Median :26.00  
##  Mean   :24.767   Mean   :26.72   Mean   :46.74   Mean   :25.92  
##  3rd Qu.:28.505   3rd Qu.:27.11   3rd Qu.:50.50   3rd Qu.:27.00  
##  Max.   :56.066   Max.   :29.20   Max.   :56.16   Max.   :30.00  
##  NA's   :537      NA's   :536     NA's   :536     NA's   :536    
##   msmt_temp_K      e_s_kPa_m       e_a_kPa_m       msmt_VPD_kPa  
##  Min.   :297.9   Min.   :3.219   Min.   :0.9894   Min.   :1.486  
##  1st Qu.:299.4   1st Qu.:3.504   1st Qu.:1.6464   1st Qu.:1.767  
##  Median :299.9   Median :3.620   Median :1.7411   Median :1.853  
##  Mean   :299.9   Mean   :3.620   Mean   :1.6833   Mean   :1.937  
##  3rd Qu.:300.3   3rd Qu.:3.701   3rd Qu.:1.7992   3rd Qu.:2.012  
##  Max.   :302.4   Max.   :4.194   Max.   :1.9326   Max.   :3.021  
##  NA's   :536     NA's   :536     NA's   :536      NA's   :536    
##       SMI         temp_mean_tmttrial temp_SD_tmttrial humidity_mean_tmttrial
##  Min.   : 7.317   Min.   :23.30      Min.   :0.5966   Min.   :13.75         
##  1st Qu.: 9.748   1st Qu.:24.05      1st Qu.:0.7828   1st Qu.:29.21         
##  Median :10.624   Median :29.74      Median :1.0461   Median :45.24         
##  Mean   :10.607   Mean   :29.61      Mean   :1.1502   Mean   :52.95         
##  3rd Qu.:11.348   3rd Qu.:35.05      3rd Qu.:1.5191   3rd Qu.:82.84         
##  Max.   :14.263   Max.   :36.00      Max.   :1.8447   Max.   :93.15         
##                                                                             
##  humidity_SD_tmttrial e_a_mean_tmttrial e_a_SD_tmttrial  VPD_mean_tmttrial
##  Min.   : 4.370       Min.   :0.4258    Min.   :0.1891   Min.   :0.2053   
##  1st Qu.: 6.234       1st Qu.:1.3896    1st Qu.:0.2651   1st Qu.:0.8018   
##  Median : 7.382       Median :2.2917    Median :0.3166   Median :2.0450   
##  Mean   : 8.758       Mean   :2.3359    Mean   :0.3934   Mean   :2.0051   
##  3rd Qu.:12.297       3rd Qu.:2.6749    3rd Qu.:0.5019   3rd Qu.:3.1651   
##  Max.   :19.846       Max.   :4.8058    Max.   :0.7691   Max.   :4.0685   
##                                                                           
##  VPD_SD_tmttrial      temp_C      humidity_percent    e_a_kPa     
##  Min.   :0.1690   Min.   :23.80   Min.   :17.90    Min.   :0.500  
##  1st Qu.:0.2558   1st Qu.:23.80   1st Qu.:34.10    1st Qu.:2.000  
##  Median :0.3483   Median :29.65   Median :56.20    Median :2.150  
##  Mean   :0.4036   Mean   :29.65   Mean   :52.85    Mean   :2.322  
##  3rd Qu.:0.4735   3rd Qu.:35.50   3rd Qu.:78.30    3rd Qu.:2.300  
##  Max.   :0.8863   Max.   :35.50   Max.   :80.90    Max.   :4.500  
##                                                                   
##     VPD_kPa     
##  Min.   :0.600  
##  1st Qu.:0.600  
##  Median :1.800  
##  Mean   :2.003  
##  3rd Qu.:3.800  
##  Max.   :3.800  
## 
dat_no_rehab %>%
  group_by(type) %>%
  summarise(mean(msmt_temp_C, na.rm = T),
            sd(msmt_temp_C, na.rm = T),
            mean(msmt_RH_percent, na.rm = T),
            sd(msmt_RH_percent, na.rm = T),
            mean(msmt_VPD_kPa, na.rm = T),
            sd(msmt_VPD_kPa, na.rm = T))
## # A tibble: 1 × 7
##   type  `mean(msmt_temp_C, na.rm = T)` sd(msmt…¹ mean(…² sd(ms…³ mean(…⁴ sd(ms…⁵
##   <fct>                          <dbl>     <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 exp                             26.7     0.799    46.7    6.76    1.94   0.337
## # … with abbreviated variable names ¹​`sd(msmt_temp_C, na.rm = T)`,
## #   ²​`mean(msmt_RH_percent, na.rm = T)`, ³​`sd(msmt_RH_percent, na.rm = T)`,
## #   ⁴​`mean(msmt_VPD_kPa, na.rm = T)`, ⁵​`sd(msmt_VPD_kPa, na.rm = T)`

Means by Day

Calculate mean values per day per tmt group.

means <- dat_no_rehab %>%
  group_by(day_n, tmt) %>%
  summarise(n_lizards = n(),
            mean_CEWL = mean(CEWL_g_m2h_mean, na.rm = T),
            sd_CEWL = sd(CEWL_g_m2h_mean, na.rm = T),
            mean_osml = mean(osmolality_mmol_kg_mean, na.rm = T),
            sd_osml = sd(osmolality_mmol_kg_mean, na.rm = T),
            mean_hct = mean(hematocrit_percent, na.rm = T),
            sd_hct = sd(hematocrit_percent, na.rm = T),
            mean_mass = mean(mass_g, na.rm = T),
            sd_mass = sd(mass_g, na.rm = T),
            mean_SMI = mean(SMI, na.rm = T),
            sd_SMI = sd(SMI, na.rm = T)) %>%
  mutate(se_CEWL = (sd_CEWL/sqrt(n_lizards)),
         se_osml = (sd_osml/sqrt(n_lizards)),
         se_hct = (sd_hct/sqrt(n_lizards)),
         se_mass = (sd_mass/sqrt(n_lizards)),
         se_SMI = (sd_SMI/sqrt(n_lizards)))
## `summarise()` has grouped output by 'day_n'. You can override using the
## `.groups` argument.
# get rid of non-defined points
# these are from when not all variables were measured for a given date
means$mean_CEWL[is.nan(means$mean_CEWL)] <- NA
means$mean_osml[is.nan(means$mean_osml)] <- NA
means$mean_hct[is.nan(means$mean_hct)] <- NA
means$mean_mass[is.nan(means$mean_mass)] <- NA
means$mean_SMI[is.nan(means$mean_SMI)] <- NA
means
## # A tibble: 24 × 18
## # Groups:   day_n [6]
##    day_n tmt      n_liz…¹ mean_…² sd_CEWL mean_…³ sd_osml mean_…⁴ sd_hct mean_…⁵
##    <dbl> <fct>      <int>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>   <dbl>
##  1     0 Cool Hu…      34    20.9    4.78    351.    20.3    39.6   5.30    11.6
##  2     0 Hot Hum…      33    21.4    4.85    347.    18.7    37.9   5.46    11.6
##  3     0 Cool Dry      33    20.0    6.08    346.    20.6    39.3   5.96    11.8
##  4     0 Hot Dry       34    20.9    5.93    347.    17.9    38.9   5.05    12.0
##  5     4 Cool Hu…      34    NA     NA       356.    24.5    34.5   5.34    11.1
##  6     4 Hot Hum…      33    NA     NA       359.    22.5    32.0   5.38    10.6
##  7     4 Cool Dry      33    NA     NA       355.    27.0    35.0   7.02    11.1
##  8     4 Hot Dry       34    NA     NA       361.    25.8    33.1   5.00    10.6
##  9     5 Cool Hu…      34    NA     NA        NA     NA      NA    NA       10.8
## 10     5 Hot Hum…      33    NA     NA        NA     NA      NA    NA       10.2
## # … with 14 more rows, 8 more variables: sd_mass <dbl>, mean_SMI <dbl>,
## #   sd_SMI <dbl>, se_CEWL <dbl>, se_osml <dbl>, se_hct <dbl>, se_mass <dbl>,
## #   se_SMI <dbl>, and abbreviated variable names ¹​n_lizards, ²​mean_CEWL,
## #   ³​mean_osml, ⁴​mean_hct, ⁵​mean_mass
# get only means for the very end
end_means <- means %>%
  dplyr::filter(day_n == 8)
write.csv(end_means, "./results_statistics/exp_end_means.csv")

End Values Only

Select all raw measurements for only day=8 values.

end_vals <- dat %>%
  dplyr::filter(day_n == 8)

delta CEWL

Get a df that only has complete observations that include CEWL values (only obs from before and after the experiment). Then, calculate the change (delta) in CEWL from before to after the experiment. Because we only measured CEWL at those two time points, it makes more sense to assess the amount of change in CEWL for each lizard, rather than measuring the change over time.

start_CEWL <- dat_no_rehab %>%
  dplyr::filter(day_n == 0) %>%
  dplyr::select(individual_ID, start_CEWL = CEWL_g_m2h_mean)
dat_no_rehab_deltaCEWL <- dat_no_rehab %>% # initiate new df
  dplyr::filter(complete.cases(CEWL_g_m2h_mean)) %>% # only use obs incl CEWL
  dplyr::filter(day_n == 8) %>% # get only obs for post-exp
  left_join(start_CEWL, by = 'individual_ID') %>% # add start CEWL to both obs for each lizard
  mutate(delta_CEWL = CEWL_g_m2h_mean - start_CEWL) # calculate deltaCEWL after-before experiment

Experiment Models

We predicted that there would be effects of day, humidity treatment, temperature treatment, and treatment VPD. However, we can’t use the standard backwards model selection because the three treatment variables are collinear (VIF much higher than acceptable) and it leads to issues with changing the sign of the estimates when all three are included together. So, we will run singular models with each treatment variable alone:

response ~ dayhumidity response ~ daytemperature response ~ day*VPD

I added models with water vapor pressure based on comments received during review, but deemed it uninformative to add to the results presented in the paper

Then, we will use AIC, RMSE, and R-sq to assess which treatment effect is most important to that response variable.

Body Mass

Building

Build each treatment effect model.

mass_mod_WVP <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * e_a_kPa +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')
mass_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')
## Warning: Model failed to converge with 1 negative eigenvalue: -1.7e+01
mass_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')
mass_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')

Assumptions

Check linear regression assumptions/conditions.

# distribution of SMI
hist(dat_no_rehab$mass_g)

# VPD model
plot(mass_mod_VPD)

simple.eda(residuals(mass_mod_VPD))

shapiro.test(residuals(mass_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mass_mod_VPD)
## W = 0.93331, p-value < 2.2e-16
# humidity model
plot(mass_mod_hum)

simple.eda(residuals(mass_mod_hum))

shapiro.test(residuals(mass_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mass_mod_hum)
## W = 0.93039, p-value < 2.2e-16
# temperature model
plot(mass_mod_temp)

simple.eda(residuals(mass_mod_temp))

shapiro.test(residuals(mass_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mass_mod_temp)
## W = 0.90802, p-value < 2.2e-16

Normality is violated, but linearity, equal error variance, and independence are all good.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

We…

  • calculate RMSE manually
  • use the r.squaredGLMM function in the MuMIn package to get the marginal R^2, which is how much of the total variance is explained by fixed effects
  • use the aictab function in the AICmodavg package to get AIC and deltaAIC values
  • get the sum of squares, F, and p-values for each variable from the anova table for each model
# calculate RMSE & R^2
mass_RMSE_Rsq <- data.frame(model = 
                               c('Day * WVP',
                                 'Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(mass_mod_WVP))^2)),
                                    sqrt(mean((residuals(mass_mod_VPD))^2)),
                                    sqrt(mean((residuals(mass_mod_hum))^2)),
                                    sqrt(mean((residuals(mass_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(mass_mod_WVP)[,"R2m"],
                                   MuMIn::r.squaredGLMM(mass_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(mass_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(mass_mod_temp)[,"R2m"]))
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
# calculate AIC
mass_models <- list(mass_mod_WVP, mass_mod_VPD, mass_mod_hum, mass_mod_temp)
EXP_mod_names <- data.frame(model = 
                               c('Day * WVP',
                                 'Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ))
mass_AICc <- data.frame(aictab(cand.set = mass_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = mass_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
mass_across <- mass_RMSE_Rsq %>%
  left_join(mass_AICc, by = 'model') %>%
  mutate(response = "Body Mass (g)") %>%
  arrange(Delta_AICc)

# calculate F & p-values
mass_VPD_p <- data.frame(anova(mass_mod_VPD, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
mass_hum_p <- data.frame(anova(mass_mod_hum, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
mass_temp_p <- data.frame(anova(mass_mod_temp, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
mass_within <- mass_VPD_p %>%
  rbind(mass_hum_p) %>%
  rbind(mass_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Body Mass (g)")

Body Condition

not in publication

Building

Build each treatment effect model.

SMI_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
SMI_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
SMI_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))

Assumptions

Check linear regression assumptions/conditions.

# distribution of SMI
hist(dat_no_rehab$SMI)

# VPD model
plot(SMI_mod_VPD)

simple.eda(residuals(SMI_mod_VPD))

shapiro.test(residuals(SMI_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(SMI_mod_VPD)
## W = 0.96044, p-value = 6.428e-14
# humidity model
plot(SMI_mod_hum)

simple.eda(residuals(SMI_mod_hum))

shapiro.test(residuals(SMI_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(SMI_mod_hum)
## W = 0.95569, p-value = 7.603e-15
# temperature model
plot(SMI_mod_temp)

simple.eda(residuals(SMI_mod_temp))

shapiro.test(residuals(SMI_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(SMI_mod_temp)
## W = 0.93877, p-value < 2.2e-16

Normality is violated, but linearity, equal error variance, and independence are all good.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

We…

  • calculate RMSE manually
  • use the r.squaredGLMM function in the MuMIn package to get the marginal R^2, which is how much of the total variance is explained by fixed effects
  • use the aictab function in the AICmodavg package to get AIC and deltaAIC values
  • get the sum of squares, F, and p-values for each variable from the anova table for each model
# calculate RMSE & R^2
SMI_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(SMI_mod_VPD))^2)),
                                    sqrt(mean((residuals(SMI_mod_hum))^2)),
                                    sqrt(mean((residuals(SMI_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(SMI_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(SMI_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(SMI_mod_temp)[,"R2m"]))

# calculate AIC
SMI_models <- list(SMI_mod_VPD, SMI_mod_hum, SMI_mod_temp)
EXP_mod_names <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ))
SMI_AICc <- data.frame(aictab(cand.set = SMI_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = SMI_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
SMI_across <- SMI_RMSE_Rsq %>%
  left_join(SMI_AICc, by = 'model') %>%
  mutate(response = "Body Condition (g')") %>%
  arrange(Delta_AICc)

# calculate F & p-values
SMI_VPD_p <- data.frame(anova(SMI_mod_VPD, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
SMI_hum_p <- data.frame(anova(SMI_mod_hum, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
SMI_temp_p <- data.frame(anova(SMI_mod_temp, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
SMI_within <- SMI_VPD_p %>%
  rbind(SMI_hum_p) %>%
  rbind(SMI_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Body Condition (g')")

Hematocrit

Building

Build each treatment effect model.

hct_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                               hematocrit_percent ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
hct_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                              hematocrit_percent ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
hct_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                              hematocrit_percent ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))

Assumptions

Check linear regression assumptions/conditions.

# distribution of 
hist(dat_no_rehab$hematocrit_percent)

# VPD model
plot(hct_mod_VPD)

simple.eda(residuals(hct_mod_VPD))

shapiro.test(residuals(hct_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_mod_VPD)
## W = 0.99652, p-value = 0.5454
# humidity model
plot(hct_mod_hum)

simple.eda(residuals(hct_mod_hum))

shapiro.test(residuals(hct_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_mod_hum)
## W = 0.99619, p-value = 0.4584
# temperature model
plot(hct_mod_temp)

simple.eda(residuals(hct_mod_temp))

shapiro.test(residuals(hct_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_mod_temp)
## W = 0.99478, p-value = 0.1975

All assumptions, normality, linearity, equal error variance, and independence are all good.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

# calculate RMSE & R^2
hct_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(hct_mod_VPD))^2)),
                                    sqrt(mean((residuals(hct_mod_hum))^2)),
                                    sqrt(mean((residuals(hct_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(hct_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(hct_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(hct_mod_temp)[,"R2m"]))

# calculate AIC
hct_models <- list(hct_mod_VPD, hct_mod_hum, hct_mod_temp)
hct_AICc <- data.frame(aictab(cand.set = hct_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = hct_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
hct_across <- hct_RMSE_Rsq %>%
  left_join(hct_AICc, by = 'model') %>%
  mutate(response = "Hematocrit (%)") %>%
  arrange(Delta_AICc)

# calculate F & p-values
hct_VPD_p <- data.frame(anova(hct_mod_VPD, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
hct_hum_p <- data.frame(anova(hct_mod_hum, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
hct_temp_p <- data.frame(anova(hct_mod_temp, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
hct_within <- hct_VPD_p %>%
  rbind(hct_hum_p) %>%
  rbind(hct_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Hematocrit (%)")

Osmolality

Building

Build each treatment effect model.

osml_mod_VPD <- lmerTest::lmer(data = dat_no_rehab,
                               osmolality_mmol_kg_mean ~ day_n * VPD_kPa +
                              (1|trial_number/individual_ID))
osml_mod_hum <- lmerTest::lmer(data = dat_no_rehab,
                               osmolality_mmol_kg_mean ~ day_n * humidity_tmt +
                              (1|trial_number/individual_ID))
osml_mod_temp <- lmerTest::lmer(data = dat_no_rehab,
                               osmolality_mmol_kg_mean ~ day_n * temp_tmt +
                              (1|trial_number/individual_ID))

Assumptions

Check linear regression assumptions/conditions.

# distribution of 
hist(dat_no_rehab$osmolality_mmol_kg_mean)

# VPD model
plot(osml_mod_VPD)

simple.eda(residuals(osml_mod_VPD))

shapiro.test(residuals(osml_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(osml_mod_VPD)
## W = 0.97687, p-value = 6.755e-06
# humidity model
plot(osml_mod_hum)

simple.eda(residuals(osml_mod_hum))

shapiro.test(residuals(osml_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(osml_mod_hum)
## W = 0.97746, p-value = 8.914e-06
# temperature model
plot(osml_mod_temp)

simple.eda(residuals(osml_mod_temp))

shapiro.test(residuals(osml_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(osml_mod_temp)
## W = 0.97346, p-value = 1.44e-06

Normality is violated, but linearity, equal error variance, and independence are all okay.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

# calculate RMSE & R^2
osml_RMSE_Rsq <- data.frame(model = 
                               c('Day * VPD',
                                 'Day * Humidity',
                                 'Day * Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(osml_mod_VPD))^2)),
                                    sqrt(mean((residuals(osml_mod_hum))^2)),
                                    sqrt(mean((residuals(osml_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(osml_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(osml_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(osml_mod_temp)[,"R2m"]))

# calculate AIC
osml_models <- list(osml_mod_VPD, osml_mod_hum, osml_mod_temp)
osml_AICc <- data.frame(aictab(cand.set = osml_models, 
                                 modnames = EXP_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = osml_models, modnames = EXP_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
osml_across <- osml_RMSE_Rsq %>%
  left_join(osml_AICc, by = 'model') %>%
  mutate(response = "Plasma Osmolality (mmol/kg)") %>%
  arrange(Delta_AICc)

# calculate F & p-values
osml_VPD_p <- data.frame(anova(osml_mod_VPD, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * VPD',
         predictor = rownames(.))
osml_hum_p <- data.frame(anova(osml_mod_hum, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Humidity',
         predictor = rownames(.))
osml_temp_p <- data.frame(anova(osml_mod_temp, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Day * Temp',
         predictor = rownames(.))

# save within model values
osml_within <- osml_VPD_p %>%
  rbind(osml_hum_p) %>%
  rbind(osml_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "Plasma Osmolality (mmol/kg)")

Body Temp

I need to double check whether CEWL has a relationship with body temperature at the point of measurement.

ggplot(dat_no_rehab) + 
  geom_point(aes(x = CEWL_g_m2h_mean,
                 y = cloacal_temp_C,
                 color = day_n))
## Warning: Removed 537 rows containing missing values (`geom_point()`).

Test an lm of raw CEWL ~ body temp for day 8 measurements only.

body_temp_test_dat <- end_vals %>%
  dplyr::filter(complete.cases(CEWL_g_m2h_mean))

CEWL_body_temp_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                         CEWL_g_m2h_mean ~ cloacal_temp_C * tmt +
                           (1|trial_number))
summary(CEWL_body_temp_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ cloacal_temp_C * tmt + (1 | trial_number)
##    Data: body_temp_test_dat
## 
## REML criterion at convergence: 839.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.02426 -0.58867 -0.09161  0.46172  3.09926 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  4.977   2.231   
##  Residual                 36.926   6.077   
## Number of obs: 133, groups:  trial_number, 5
## 
## Fixed effects:
##                              Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                  42.83759   26.65061 123.21107   1.607    0.111
## cloacal_temp_C               -0.49420    1.05450 122.93003  -0.469    0.640
## tmtHot Humid                -16.46275   34.17703 122.09739  -0.482    0.631
## tmtCool Dry                 -23.23095   47.80497 123.10602  -0.486    0.628
## tmtHot Dry                   -3.95989   36.73317 123.70420  -0.108    0.914
## cloacal_temp_C:tmtHot Humid   0.90966    1.35821 122.13178   0.670    0.504
## cloacal_temp_C:tmtCool Dry    0.65396    1.86629 123.06213   0.350    0.727
## cloacal_temp_C:tmtHot Dry    -0.06271    1.44351 123.65911  -0.043    0.965
## 
## Correlation of Fixed Effects:
##             (Intr) clc__C tmtHtH tmtClD tmtHtD c__C:H c__C:CD
## clocl_tmp_C -0.999                                           
## tmtHotHumid -0.745  0.744                                    
## tmtCool Dry -0.528  0.528  0.434                             
## tmtHot Dry  -0.694  0.694  0.564  0.422                      
## clcl_t_C:HH  0.741 -0.741 -0.999 -0.432 -0.562               
## clcl_t_C:CD  0.536 -0.536 -0.439 -0.999 -0.427  0.438        
## clcl_t_C:HD  0.699 -0.700 -0.567 -0.425 -0.999  0.566  0.430
anova(CEWL_body_temp_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                     Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## cloacal_temp_C      1.1914  1.1914     1 115.35  0.0323 0.8578
## tmt                14.9005  4.9668     3 122.58  0.1345 0.9393
## cloacal_temp_C:tmt 27.5730  9.1910     3 122.56  0.2489 0.8620
ggplot(body_temp_test_dat) + 
  geom_point(aes(x = msmt_temp_C,
                 y = cloacal_temp_C)) +
  xlim(23, 28) + ylim(23,28)

ggplot(body_temp_test_dat) + 
  geom_point(aes(x = msmt_temp_C,
                 y = CEWL_g_m2h_mean,
                 color = tmt))

Do the same exact stats for ambient temp and VPD at the time of measurement:

CEWL_msmt_temp_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                                     CEWL_g_m2h_mean ~ msmt_temp_C * tmt +
                                       (1|trial_number))
CEWL_msmt_WVP_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                                     CEWL_g_m2h_mean ~ e_a_kPa_m * tmt +
                                       (1|trial_number))
CEWL_msmt_VPD_mod <- lmerTest::lmer(data = body_temp_test_dat, 
                                     CEWL_g_m2h_mean ~ msmt_VPD_kPa * tmt +
                                       (1|trial_number))

summary(CEWL_msmt_temp_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ msmt_temp_C * tmt + (1 | trial_number)
##    Data: body_temp_test_dat
## 
## REML criterion at convergence: 824.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8115 -0.5673 -0.0847  0.4720  3.5376 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  7.298   2.701   
##  Residual                 33.499   5.788   
## Number of obs: 133, groups:  trial_number, 5
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)              -39.6305    48.2297 112.4351  -0.822    0.413
## msmt_temp_C                2.6733     1.8380 112.5242   1.454    0.149
## tmtHot Humid             -73.4611    65.2521 121.4757  -1.126    0.262
## tmtCool Dry              -18.6586    60.6146 121.3298  -0.308    0.759
## tmtHot Dry               -61.0925    64.0861 121.2408  -0.953    0.342
## msmt_temp_C:tmtHot Humid   3.0314     2.4837 121.4794   1.221    0.225
## msmt_temp_C:tmtCool Dry    0.4532     2.3073 121.3153   0.196    0.845
## msmt_temp_C:tmtHot Dry     2.1152     2.4428 121.2415   0.866    0.388
## 
## Correlation of Fixed Effects:
##             (Intr) msm__C tmtHtH tmtClD tmtHtD m__C:H m__C:CD
## msmt_temp_C -0.999                                           
## tmtHotHumid -0.593  0.592                                    
## tmtCool Dry -0.614  0.613  0.467                             
## tmtHot Dry  -0.580  0.580  0.445  0.478                      
## msmt_t_C:HH  0.594 -0.594 -1.000 -0.467 -0.445               
## msmt_t_C:CD  0.614 -0.614 -0.467 -1.000 -0.478  0.467        
## msmt_t_C:HD  0.579 -0.579 -0.444 -0.478 -1.000  0.445  0.478
summary(CEWL_msmt_VPD_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ msmt_VPD_kPa * tmt + (1 | trial_number)
##    Data: body_temp_test_dat
## 
## REML criterion at convergence: 817.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8187 -0.5552 -0.0438  0.3596  3.2064 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 11.75    3.428   
##  Residual                 34.18    5.846   
## Number of obs: 133, groups:  trial_number, 5
## 
## Fixed effects:
##                           Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)                 21.583     14.417  63.780   1.497   0.1393  
## msmt_VPD_kPa                 4.952      8.054  67.896   0.615   0.5407  
## tmtHot Humid               -24.082     16.730 120.618  -1.439   0.1526  
## tmtCool Dry                 -9.284     16.225 120.721  -0.572   0.5682  
## tmtHot Dry                 -25.308     16.211 120.503  -1.561   0.1211  
## msmt_VPD_kPa:tmtHot Humid   17.094      9.373 120.640   1.824   0.0707 .
## msmt_VPD_kPa:tmtCool Dry     1.504      9.025 120.677   0.167   0.8679  
## msmt_VPD_kPa:tmtHot Dry     11.097      9.093 120.529   1.220   0.2247  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) m_VPD_ tmtHtH tmtClD tmtHtD m_VPDH m_VPD_P:CD
## msmt_VPD_kP -0.992                                              
## tmtHotHumid -0.538  0.538                                       
## tmtCool Dry -0.524  0.523  0.472                                
## tmtHot Dry  -0.535  0.535  0.476  0.490                         
## ms_VPD_P:HH  0.538 -0.542 -0.996 -0.470 -0.474                  
## ms_VPD_P:CD  0.527 -0.531 -0.474 -0.996 -0.492  0.476           
## ms_VPD_P:HD  0.529 -0.533 -0.473 -0.488 -0.996  0.475  0.493
anova(CEWL_msmt_temp_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                  Sum Sq Mean Sq NumDF   DenDF F value   Pr(>F)   
## msmt_temp_C     309.546 309.546     1  59.634  9.2405 0.003513 **
## tmt              57.484  19.161     3 121.536  0.5720 0.634497   
## msmt_temp_C:tmt  65.354  21.785     3 121.538  0.6503 0.584244   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(CEWL_msmt_WVP_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##               Sum Sq Mean Sq NumDF   DenDF F value  Pr(>F)  
## e_a_kPa_m     141.13 141.126     1   8.973  4.0770 0.07432 .
## tmt           174.69  58.230     3 121.755  1.6822 0.17438  
## e_a_kPa_m:tmt 129.98  43.325     3 121.753  1.2516 0.29414  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(CEWL_msmt_VPD_mod, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                  Sum Sq Mean Sq NumDF   DenDF F value  Pr(>F)  
## msmt_VPD_kPa     115.03 115.026     1  32.112  3.3653 0.07586 .
## tmt              112.93  37.643     3 121.388  1.1013 0.35140  
## msmt_VPD_kPa:tmt 152.90  50.966     3 121.397  1.4911 0.22039  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CEWL

First, get a separate df of end-experiment delta CEWL for each tmt group:

dat_no_rehab_deltaCEWL_post <- dat_no_rehab_deltaCEWL %>%
  dplyr::filter(day_n == 8)
HH8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Hu")
HD8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Dr")
CH8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool H")
CD8 <- dat_no_rehab_deltaCEWL_post %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool D")

t-tests

I use t-tests to assess whether the change in CEWL from the experiment is significantly different from zero for each treatment group.

CEWL_ttest_HH <- t.test(HH8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_HH
## 
##  One Sample t-test
## 
## data:  HH8$delta_CEWL
## t = 8.7341, df = 32, p-value = 5.573e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  11.76821 18.92672
## sample estimates:
## mean of x 
##  15.34746
CEWL_ttest_HD <- t.test(HD8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_HD
## 
##  One Sample t-test
## 
## data:  HD8$delta_CEWL
## t = 2.513, df = 33, p-value = 0.01703
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.7042497 6.6929758
## sample estimates:
## mean of x 
##  3.698613
CEWL_ttest_CH <- t.test(CH8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_CH
## 
##  One Sample t-test
## 
## data:  CH8$delta_CEWL
## t = 8.7488, df = 32, p-value = 5.363e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   7.296566 11.725292
## sample estimates:
## mean of x 
##  9.510929
CEWL_ttest_CD <- t.test(CD8$delta_CEWL, 
                        mu = 0, alternative = "two.sided")
CEWL_ttest_CD
## 
##  One Sample t-test
## 
## data:  CD8$delta_CEWL
## t = 2.75, df = 32, p-value = 0.009721
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.9302207 6.2446177
## sample estimates:
## mean of x 
##  3.587419

Building

Build each treatment effect model. These are equivalent to the models built for body condition, hematocrit, and plasma osmolality, but the effect of time is not included because we assess CEWL as delta CEWL.

CEWL_mod_VPD <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ VPD_kPa +
                              (1|trial_number))
CEWL_mod_WVP <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ e_a_kPa +
                              (1|trial_number))
CEWL_mod_hum <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ humidity_tmt +
                              (1|trial_number))
CEWL_mod_temp <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ temp_tmt +
                              (1|trial_number))
CEWL_mod_double <- lmerTest::lmer(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ temp_tmt*humidity_tmt +
                              (1|trial_number))

# summary(CEWL_mod_VPD)
# summary(CEWL_mod_hum)
# summary(CEWL_mod_temp)
# summary(CEWL_mod_double)

Assumptions

Check linear regression assumptions/conditions.

# distribution of 
hist(dat_no_rehab_deltaCEWL$delta_CEWL)

# VPD model
plot(CEWL_mod_VPD)

simple.eda(residuals(CEWL_mod_VPD))

shapiro.test(residuals(CEWL_mod_VPD))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_mod_VPD)
## W = 0.97248, p-value = 0.008418
# humidity model
plot(CEWL_mod_hum)

simple.eda(residuals(CEWL_mod_hum))

shapiro.test(residuals(CEWL_mod_hum))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_mod_hum)
## W = 0.98319, p-value = 0.1002
# temperature model
plot(CEWL_mod_temp)

simple.eda(residuals(CEWL_mod_temp))

shapiro.test(residuals(CEWL_mod_temp))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_mod_temp)
## W = 0.98262, p-value = 0.08761

All assumptions are fine.

Comparison

Now, compare the AIC, RMSE, and R^2 values across models, and the F and p values of the variables for each model.

# calculate RMSE & R^2
CEWL_RMSE_Rsq <- data.frame(model = 
                               c('WVP',
                                 'VPD',
                                 'Humidity',
                                 'Temp'
                                 ),
                           RMSE = c(sqrt(mean((residuals(CEWL_mod_WVP))^2)),
                                    sqrt(mean((residuals(CEWL_mod_VPD))^2)),
                                    sqrt(mean((residuals(CEWL_mod_hum))^2)),
                                    sqrt(mean((residuals(CEWL_mod_temp))^2))),
                           # marginal Rsq for the amount of variance
                           # explained by fixed effects only
                           Rsq = c(MuMIn::r.squaredGLMM(CEWL_mod_WVP)[,"R2m"],
                                   MuMIn::r.squaredGLMM(CEWL_mod_VPD)[,"R2m"],
                                   MuMIn::r.squaredGLMM(CEWL_mod_hum)[,"R2m"],
                                   MuMIn::r.squaredGLMM(CEWL_mod_temp)[,"R2m"]))

# calculate AIC
CEWL_models <- list(CEWL_mod_WVP, CEWL_mod_VPD, CEWL_mod_hum, CEWL_mod_temp)
CEWL_mod_names <- data.frame(model = 
                               c('WVP',
                                 'VPD',
                                 'Humidity',
                                 'Temp'
                                 ))
CEWL_AICc <- data.frame(aictab(cand.set = CEWL_models, 
                                 modnames = CEWL_mod_names))
## Warning in aictab.AIClmerModLmerTest(cand.set = CEWL_models, modnames = CEWL_mod_names): 
## Model selection for fixed effects is only appropriate with ML estimation:
## REML (default) should only be used to select random effects for a constant set of fixed effects
# compare across models
CEWL_across <- CEWL_RMSE_Rsq %>%
  left_join(CEWL_AICc, by = 'model') %>%
  mutate(response = "deltaCEWL") %>%
  arrange(Delta_AICc)

# calculate F & p-values
CEWL_VPD_p <- data.frame(anova(CEWL_mod_VPD, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'VPD',
         predictor = rownames(.))
CEWL_hum_p <- data.frame(anova(CEWL_mod_hum, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Humidity',
         predictor = rownames(.))
CEWL_temp_p <- data.frame(anova(CEWL_mod_temp, 
                              type = "3", 
                              ddf = "Kenward-Roger")) %>%
  mutate(model = 'Temp',
         predictor = rownames(.))

# save within model values
CEWL_within <- CEWL_VPD_p %>%
  rbind(CEWL_hum_p) %>%
  rbind(CEWL_temp_p) %>%
  mutate(df = paste((NumDF), round(DenDF, 0), sep = ", "),
         response = "deltaCEWL")

But, water vapor pressure is pretty confounded with humidity tmt.

Group Export

Put all the model statistics into one df/csv - one for among-model comparisons, and one for within-model stats.

# comparisons of all models for each variable
experiment_model_compare <- CEWL_across %>%
  rbind(osml_across) %>%
  rbind(hct_across) %>%
  rbind(mass_across) %>%
  rbind(SMI_across) %>%
  dplyr::select(response, model, 
                RMSE, Rsq, AICc, Delta_AICc) %>%
  mutate(RMSE = round(RMSE, 2), 
         Rsq = round(Rsq, 2), 
         AICc = round(AICc, 2), 
         Delta_AICc = round(Delta_AICc, 2)) 
write.csv(experiment_model_compare, 
          "./results_statistics/exp_model_comparisons.csv")

# all of the actual model results
experiment_model_values <- CEWL_within %>%
  rbind(osml_within) %>%
  rbind(hct_within) %>%
  rbind(mass_within) %>%
  rbind(SMI_within) %>%
  dplyr::select(response, model, predictor, 
                seq_sum_of_squares = Sum.Sq, 
                df, 
                F_statistic = F.value, 
                p_value = Pr..F.) %>%
  mutate(seq_sum_of_squares = round(seq_sum_of_squares, 0), 
         F_statistic = round(F_statistic, 2)) 
write.csv(experiment_model_values, "./results_statistics/exp_model_values.csv")

Effect Estimates

End Value CIs

Now, we can use the emmeans function from the emmeans package to estimate marginal means and confidence intervals for the values among treatment groups at the end of the experiment. But, to get for each treatment group, we need to run a new model with day 8 data only and tmt as a single, 4-category variable.

# Body Mass
mass_mod_end <- lmerTest::lmer(data = end_vals,
                              mass_g ~ tmt + 
                              (1|trial_number))
## boundary (singular) fit: see help('isSingular')
mass_emmeans <- data.frame(emmeans(mass_mod_end, "tmt")) %>%
  mutate(response = "Body Mass (g)")
mass_pairwise <- data.frame(pairs(emmeans(mass_mod_end, "tmt"))) %>%
  mutate(response = "Body Mass (g)")

# Body Condition
SMI_mod_end <- lmerTest::lmer(data = end_vals,
                              SMI ~ tmt + 
                              (1|trial_number))
SMI_emmeans <- data.frame(emmeans(SMI_mod_end, "tmt")) %>%
  mutate(response = "Body Condition (g')")
SMI_pairwise <- data.frame(pairs(emmeans(SMI_mod_end, "tmt"))) %>%
  mutate(response = "Body Condition (g')")

# Hematocrit
hct_mod_end <- lmerTest::lmer(data = end_vals,
                              hematocrit_percent ~ tmt + 
                              (1|trial_number))
hct_emmeans <- data.frame(emmeans(hct_mod_end, "tmt")) %>%
  mutate(response = "Hematocrit (%)")
hct_pairwise <- data.frame(pairs(emmeans(hct_mod_end, "tmt"))) %>%
  mutate(response = "Hematocrit (%)")

# Plasma Osmolality
osml_mod_end <- lmerTest::lmer(data = end_vals,
                              osmolality_mmol_kg_mean ~ tmt + 
                              (1|trial_number))
osml_emmeans <- data.frame(emmeans(osml_mod_end, "tmt")) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")
osml_pairwise <- data.frame(pairs(emmeans(osml_mod_end, "tmt"))) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")

# CEWL
CEWL_mod_end <- lmerTest::lmer(data = end_vals,
                              CEWL_g_m2h_mean ~ tmt + 
                              (1|trial_number))
CEWL_emmeans <- data.frame(emmeans(CEWL_mod_end, "tmt")) %>%
  mutate(response = "CEWL (g/m2h)")
CEWL_pairwise <- data.frame(pairs(emmeans(CEWL_mod_end, "tmt"))) %>%
  mutate(response = "CEWL (g/m2h)")

Rate of Change

# Body MASS
mass_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              mass_g ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
## boundary (singular) fit: see help('isSingular')
mass_emtrends <- data.frame(emtrends(mass_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Body Mass (g)")
mass_pairtrend <- data.frame(pairs(emtrends(mass_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Body Mass (g)")

# Body Condition
SMI_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              SMI ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
SMI_emtrends <- data.frame(emtrends(SMI_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Body Condition (g')")
SMI_pairtrend <- data.frame(pairs(emtrends(SMI_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Body Condition (g')")

# Hematocrit
hct_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              hematocrit_percent ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
hct_emtrends <- data.frame(emtrends(hct_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Hematocrit (%)")
hct_pairtrend <- data.frame(pairs(emtrends(hct_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Hematocrit (%)")

# Plasma Osmolality
osml_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              osmolality_mmol_kg_mean ~ day_n * tmt + 
                              (1|trial_number/individual_ID))
osml_emtrends <- data.frame(emtrends(osml_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")
osml_pairtrend <- data.frame(pairs(emtrends(osml_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "Plasma Osmolality (mmol/kg)")

# CEWL
CEWL_mod_day <- lmerTest::lmer(data = dat_no_rehab,
                              CEWL_g_m2h_mean ~ day_n * tmt + 
                              (1|trial_number))
CEWL_emtrends <- data.frame(emtrends(CEWL_mod_day, "tmt", var = "day_n")) %>%
  mutate(response = "CEWL (g/m2h)")
CEWL_pairtrend <- data.frame(pairs(emtrends(CEWL_mod_day, "tmt", var = "day_n"))) %>%
  mutate(response = "CEWL (g/m2h)")

# put together
all_emtrends <- CEWL_emtrends %>%
  rbind(osml_emtrends) %>%
  rbind(hct_emtrends) %>%
  rbind(mass_emtrends) %>%
  rbind(SMI_emtrends) %>%
  mutate(confint95 = paste(round(lower.CL, 2), round(upper.CL, 2), sep = ", ")) %>%
  dplyr::select(response, tmt, 
                per_day_change = day_n.trend,
                confint95,
                SE, df)
write.csv(all_emtrends, "./results_statistics/exp_emtrends_per_day_change.csv")

Given the daily trends for plasma osmolality, in total for acclimation change in osml had 1.04 (0.13x8) to 13 (1.587x8) mmol/kg of change.

CEWL ~ osmolality

Let’s compare the relationship of CEWL ~ osmolality post-experiment to what we found pre-experiment in the “analysis_capture” file. The treatment-specific sub-dataframes were created in the “Body Temp” model subsection above.

These are NOT the stats presents for Figure 6

CEWL_osml <- lmerTest::lmer(data = end_vals, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: end_vals
## 
## REML criterion at convergence: 855.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1565 -0.6405 -0.2787  0.6163  3.5785 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  5.66    2.379   
##  Residual                 58.25    7.632   
## Number of obs: 123, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)             18.53375   10.95312 66.51336   1.692   0.0953 .
## osmolality_mmol_kg_mean  0.02904    0.03081 72.84361   0.942   0.3491  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.993
anova(CEWL_osml, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 45.195  45.195     1 69.109  0.7759 0.3814
HH_CEWL_osml <- lmerTest::lmer(data = HH8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
## boundary (singular) fit: see help('isSingular')
summary(HH_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: HH8
## 
## REML criterion at convergence: 222.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1357 -0.6398 -0.1682  0.4261  2.4440 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  0.00    0.00    
##  Residual                 62.57    7.91    
## Number of obs: 32, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)             19.77938   23.68747 30.00000   0.835    0.410
## osmolality_mmol_kg_mean  0.04658    0.06734 30.00000   0.692    0.494
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.998
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
anova(HH_CEWL_osml, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 24.554  24.554     1 18.124  0.3924 0.5388
HD_CEWL_osml <- lmerTest::lmer(data = HD8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(HD_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: HD8
## 
## REML criterion at convergence: 171.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55684 -0.59339 -0.07624  0.40035  2.96626 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  6.99    2.644   
##  Residual                 11.28    3.358   
## Number of obs: 31, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)              2.11489   11.95116 25.69290   0.177    0.861  
## osmolality_mmol_kg_mean  0.06101    0.03291 26.24328   1.854    0.075 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.994
anova(HD_CEWL_osml, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 31.917  31.917     1 26.042  2.8298 0.1045
CH_CEWL_osml <- lmerTest::lmer(data = CH8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(CH_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: CH8
## 
## REML criterion at convergence: 192.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2869 -0.7878 -0.1247  0.6553  2.0089 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept)  6.866   2.620   
##  Residual                 30.995   5.567   
## Number of obs: 30, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)             22.02423   13.70141 27.96594   1.607    0.119
## osmolality_mmol_kg_mean  0.02489    0.03901 27.97274   0.638    0.529
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.993
anova(CH_CEWL_osml, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## osmolality_mmol_kg_mean 11.376  11.376     1 27.97   0.367 0.5495
CD_CEWL_osml <- lmerTest::lmer(data = CD8, 
                               CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean +
                               (1|trial_number))
summary(CD_CEWL_osml)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h_mean ~ osmolality_mmol_kg_mean + (1 | trial_number)
##    Data: CD8
## 
## REML criterion at convergence: 170
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8562 -0.5019  0.1381  0.4084  1.8415 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  trial_number (Intercept) 21.14    4.598   
##  Residual                 11.34    3.368   
## Number of obs: 30, groups:  trial_number, 5
## 
## Fixed effects:
##                         Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)             -1.45589   11.16955 27.92145  -0.130   0.8972  
## osmolality_mmol_kg_mean  0.07357    0.03133 27.93722   2.348   0.0262 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## osmllty_m__ -0.981
anova(CD_CEWL_osml, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##                         Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## osmolality_mmol_kg_mean 55.393  55.393     1 27.937  4.8837 0.03548 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Even though only data for one treatment group was significant, all have a positive relationship

CEWL ~ WVP

The Rsq is better… 0.22 vs 0.15

CEWL_WVP_lm <- lm(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ e_a_mean_tmttrial)
plot(CEWL_WVP_lm)

simple.eda(residuals(CEWL_WVP_lm))

shapiro.test(residuals(CEWL_WVP_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_WVP_lm)
## W = 0.98635, p-value = 0.209
summary(CEWL_WVP_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ e_a_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.6389  -6.4501  -0.5578   4.6940  22.9179 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.7748     1.3846   0.560    0.577    
## e_a_mean_tmttrial   3.0979     0.5049   6.136 9.37e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.388 on 131 degrees of freedom
## Multiple R-squared:  0.2232, Adjusted R-squared:  0.2173 
## F-statistic: 37.65 on 1 and 131 DF,  p-value: 9.373e-09
CEWL_WVP_poly <- lm(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ poly(e_a_mean_tmttrial, 2))
plot(CEWL_WVP_poly)

simple.eda(residuals(CEWL_WVP_poly))

shapiro.test(residuals(CEWL_WVP_poly))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_WVP_poly)
## W = 0.98651, p-value = 0.2166
summary(CEWL_WVP_poly)
## 
## Call:
## lm(formula = delta_CEWL ~ poly(e_a_mean_tmttrial, 2), data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.3251  -6.0525  -0.4557   5.0184  22.4236 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   8.0035     0.7271  11.008  < 2e-16 ***
## poly(e_a_mean_tmttrial, 2)1  51.4627     8.3851   6.137 9.44e-09 ***
## poly(e_a_mean_tmttrial, 2)2   8.7057     8.3851   1.038    0.301    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.385 on 130 degrees of freedom
## Multiple R-squared:  0.2296, Adjusted R-squared:  0.2178 
## F-statistic: 19.37 on 2 and 130 DF,  p-value: 4.327e-08
anova(CEWL_WVP_poly, CEWL_WVP_lm)
## Analysis of Variance Table
## 
## Model 1: delta_CEWL ~ poly(e_a_mean_tmttrial, 2)
## Model 2: delta_CEWL ~ e_a_mean_tmttrial
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    130 9140.3                           
## 2    131 9216.1 -1   -75.788 1.0779 0.3011

CEWL ~ VPD

We also want to assess whether the vapor pressure deficit lizards experienced during acclimation has an effect on the change in CEWL.

CEWL_VPD_lm <- lm(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ VPD_mean_tmttrial)
plot(CEWL_VPD_lm)

simple.eda(residuals(CEWL_VPD_lm))

shapiro.test(residuals(CEWL_VPD_lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_VPD_lm)
## W = 0.96688, p-value = 0.002494
summary(CEWL_VPD_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ VPD_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8041  -5.6085  -0.7356   5.5386  27.3599 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        13.9609     1.4186   9.841  < 2e-16 ***
## VPD_mean_tmttrial  -2.9512     0.5943  -4.965 2.09e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.731 on 131 degrees of freedom
## Multiple R-squared:  0.1584, Adjusted R-squared:  0.152 
## F-statistic: 24.66 on 1 and 131 DF,  p-value: 2.095e-06

Even though the data is slightly nonlinear, a linear model does a fine job explaining the data.

To be extra sure, also check for high leverage and influential points using H-values and Cook’s distance:

# compute values for observations 
high_leverage <- data.frame(H = hatvalues(CEWL_VPD_lm)) %>% 
  mutate(row = row_number())

# compute cutoff value 
h_bar <- (3*sum(high_leverage$H))/nrow(high_leverage)

# add to original dataframe 
# see which observations have extremely high leverage (if any)
high_leverage_dat <- dat_no_rehab_deltaCEWL %>%
  mutate(row = row_number()) %>%
  left_join(., high_leverage, by = "row") %>%
  dplyr::filter(H > h_bar) 
high_leverage_dat
## # A tibble: 0 × 38
## # Groups:   individual_ID [0]
## # … with 38 variables: measurement_date <date>, type <fct>,
## #   individual_ID <fct>, mass_g <dbl>, hematocrit_percent <int>,
## #   trial_number <fct>, temp_tmt <fct>, humidity_tmt <fct>, SVL_mm <int>,
## #   tmt <fct>, day_n <dbl>, day_factor <fct>, osmolality_mmol_kg_mean <dbl>,
## #   CEWL_g_m2h_mean <dbl>, msmt_temp_C <dbl>, msmt_RH_percent <dbl>,
## #   cloacal_temp_C <dbl>, msmt_temp_K <dbl>, e_s_kPa_m <dbl>, e_a_kPa_m <dbl>,
## #   msmt_VPD_kPa <dbl>, SMI <dbl>, temp_mean_tmttrial <dbl>, …
# get Cook's distance 
cooks <- data.frame(c = cooks.distance(CEWL_VPD_lm)) %>% 
  mutate(row = row_number())

# add to original dataframe 
influential <- dat_no_rehab_deltaCEWL %>%
  mutate(row = row_number()) %>% 
  left_join(., cooks, by = "row")

# see moderately influential points 
cook_mod_inf <- influential %>% 
  dplyr::filter(c>0.5) 
cook_mod_inf
## # A tibble: 0 × 38
## # Groups:   individual_ID [0]
## # … with 38 variables: measurement_date <date>, type <fct>,
## #   individual_ID <fct>, mass_g <dbl>, hematocrit_percent <int>,
## #   trial_number <fct>, temp_tmt <fct>, humidity_tmt <fct>, SVL_mm <int>,
## #   tmt <fct>, day_n <dbl>, day_factor <fct>, osmolality_mmol_kg_mean <dbl>,
## #   CEWL_g_m2h_mean <dbl>, msmt_temp_C <dbl>, msmt_RH_percent <dbl>,
## #   cloacal_temp_C <dbl>, msmt_temp_K <dbl>, e_s_kPa_m <dbl>, e_a_kPa_m <dbl>,
## #   msmt_VPD_kPa <dbl>, SMI <dbl>, temp_mean_tmttrial <dbl>, …

NO high leverage or influential points! :)

We will double check a comparison of a polynomial model, just to be sure:

CEWL_VPD_poly <- lm(data = dat_no_rehab_deltaCEWL,
                               delta_CEWL ~ poly(VPD_mean_tmttrial, 2))
plot(CEWL_VPD_poly)

simple.eda(residuals(CEWL_VPD_poly))

shapiro.test(residuals(CEWL_VPD_poly))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(CEWL_VPD_poly)
## W = 0.97275, p-value = 0.008935
summary(CEWL_VPD_poly)
## 
## Call:
## lm(formula = delta_CEWL ~ poly(VPD_mean_tmttrial, 2), data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.2213  -5.2906  -0.9053   4.9843  26.5733 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   8.0035     0.7523  10.639  < 2e-16 ***
## poly(VPD_mean_tmttrial, 2)1 -43.3514     8.6760  -4.997 1.84e-06 ***
## poly(VPD_mean_tmttrial, 2)2 -14.1269     8.6760  -1.628    0.106    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.676 on 130 degrees of freedom
## Multiple R-squared:  0.1752, Adjusted R-squared:  0.1625 
## F-statistic: 13.81 on 2 and 130 DF,  p-value: 3.647e-06

LINE assumptions are equally-well-met.

The polynomial factor is not significant, but the R-sq is slightly higher for the poly model.

Compare WVP / VPD / poly

Compare RMSE and AIC:

sqrt(mean((residuals(CEWL_VPD_lm))^2))
## [1] 8.664654
sqrt(mean((residuals(CEWL_VPD_poly))^2))
## [1] 8.577628
sqrt(mean((residuals(CEWL_WVP_lm))^2))
## [1] 8.324289
sqrt(mean((residuals(CEWL_WVP_poly))^2))
## [1] 8.289991
CEWL_VPD_models <- list(CEWL_WVP_lm, CEWL_WVP_poly,
                        CEWL_VPD_lm, CEWL_VPD_poly)
CEWL_VPD_mod_names <- data.frame(model = 
                               c('WVP linear',
                                 'WVP polynomial',
                                 'VPD linear',
                                 'VPD polynomial'
                                 ))
CEWL_VPD_AICc <- data.frame(aictab(cand.set = CEWL_VPD_models, 
                                 modnames = CEWL_VPD_mod_names))
CEWL_VPD_AICc
##            model K     AICc Delta_AICc    ModelLik      AICcWt        LL
## 1     WVP linear 3 947.3249   0.000000 1.000000000 0.621392996 -470.5695
## 2 WVP polynomial 4 948.3532   1.028205 0.598037054 0.371616037 -470.0203
## 4 VPD polynomial 4 957.4260  10.101091 0.006405839 0.003980544 -474.5568
## 3     VPD linear 3 957.9847  10.659766 0.004844638 0.003010424 -475.8993
##      Cum.Wt
## 1 0.6213930
## 2 0.9930090
## 4 0.9969896
## 3 1.0000000

RMSE is slightly lower for the polynomial model. But, AIC is not meaningfully different between the two versions for each variable.

I’ll use the lm for each as the best models.

summary(CEWL_VPD_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ VPD_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8041  -5.6085  -0.7356   5.5386  27.3599 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        13.9609     1.4186   9.841  < 2e-16 ***
## VPD_mean_tmttrial  -2.9512     0.5943  -4.965 2.09e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.731 on 131 degrees of freedom
## Multiple R-squared:  0.1584, Adjusted R-squared:  0.152 
## F-statistic: 24.66 on 1 and 131 DF,  p-value: 2.095e-06
anova(CEWL_VPD_lm) # single factor slr = can't specify ddf and type SS the same way (& don't need to I think)
## Analysis of Variance Table
## 
## Response: delta_CEWL
##                    Df Sum Sq Mean Sq F value    Pr(>F)    
## VPD_mean_tmttrial   1 1879.3 1879.34  24.656 2.095e-06 ***
## Residuals         131 9985.1   76.22                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(CEWL_WVP_lm)
## 
## Call:
## lm(formula = delta_CEWL ~ e_a_mean_tmttrial, data = dat_no_rehab_deltaCEWL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.6389  -6.4501  -0.5578   4.6940  22.9179 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.7748     1.3846   0.560    0.577    
## e_a_mean_tmttrial   3.0979     0.5049   6.136 9.37e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.388 on 131 degrees of freedom
## Multiple R-squared:  0.2232, Adjusted R-squared:  0.2173 
## F-statistic: 37.65 on 1 and 131 DF,  p-value: 9.373e-09
anova(CEWL_WVP_lm)
## Analysis of Variance Table
## 
## Response: delta_CEWL
##                    Df Sum Sq Mean Sq F value    Pr(>F)    
## e_a_mean_tmttrial   1 2648.4 2648.41  37.645 9.373e-09 ***
## Residuals         131 9216.1   70.35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

stats for paper

VPD: (estimate = -2.95, SE = 0.59, SS = 1879, F(1,131) = 24.66, p < 0.0001, adj-Rsq = 0.15)

WVP: (estimate = 3.10, SE = 0.50, SS = 2648, F(1,131) = 37.65, p < 0.0001, adj-Rsq = 0.22)

Recovery Models

not in publication

I want to know how the 2-day recovery period affects physiology relative to post- and pre- experiment. To do this, I’ll use two-sided t-tests comparing recovery delta values to the hypothesis of mu=0.

SMI

SMI_rmod_post_exp <- t.test(recovery_v_post_exp$delta_SMI_10_8, 
                           mu = 0, alternative = "two.sided")
SMI_rmod_post_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_post_exp$delta_SMI_10_8
## t = 6.677, df = 131, p-value = 6.292e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3171205 0.5841444
## sample estimates:
## mean of x 
## 0.4506324
SMI_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_SMI_10_0, 
                           mu = 0, alternative = "two.sided")
SMI_rmod_pre_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_pre_exp$delta_SMI_10_0
## t = -11.542, df = 131, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.3426322 -0.9497445
## sample estimates:
## mean of x 
## -1.146188

Hematocrit

hct_rmod_post_exp <- t.test(recovery_v_post_exp$delta_hct_10_8, 
                           mu = 0, alternative = "two.sided")
hct_rmod_post_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_post_exp$delta_hct_10_8
## t = -4.2083, df = 126, p-value = 4.85e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -3.033119 -1.092866
## sample estimates:
## mean of x 
## -2.062992
hct_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_hct_10_0, 
                           mu = 0, alternative = "two.sided")
hct_rmod_pre_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_pre_exp$delta_hct_10_0
## t = -22.249, df = 129, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -13.50271 -11.29729
## sample estimates:
## mean of x 
##     -12.4

Osmolality

osml_rmod_post_exp <- t.test(recovery_v_post_exp$delta_osml_10_8, 
                           mu = 0, alternative = "two.sided")
osml_rmod_post_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_post_exp$delta_osml_10_8
## t = 3.4782, df = 121, p-value = 0.0007021
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   4.330203 15.772256
## sample estimates:
## mean of x 
##  10.05123
osml_rmod_pre_exp <- t.test(recovery_v_pre_exp$delta_osml_10_0, 
                           mu = 0, alternative = "two.sided")
osml_rmod_pre_exp
## 
##  One Sample t-test
## 
## data:  recovery_v_pre_exp$delta_osml_10_0
## t = 3.75, df = 130, p-value = 0.0002651
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   5.758213 18.618378
## sample estimates:
## mean of x 
##   12.1883

Group Export

Now, I put all the recovery t-test statistics into one dataframe, and export it.

recovery_stats <- broom.mixed::tidy(osml_rmod_pre_exp) %>%
  rbind(broom.mixed::tidy(osml_rmod_post_exp)) %>%
  rbind(broom.mixed::tidy(hct_rmod_pre_exp)) %>%
  rbind(broom.mixed::tidy(hct_rmod_post_exp)) %>%
  rbind(broom.mixed::tidy(SMI_rmod_pre_exp)) %>%
  rbind(broom.mixed::tidy(SMI_rmod_post_exp)) %>%
  mutate(response = c(rep("Plasma Osmolality (mmol/kg)", 2),
                      rep("Hematocrit (%)", 2),
                      rep("Body Condition (g')", 2)),
         pre_post_exp = c(rep(c("pre", "post"), 3)))
write.csv(recovery_stats, 
          "./results_statistics/recovery_stats.csv")

Figures

Colors & Shapes

Create custom colors, shapes, and labels to make it easy to apply changes across all figures.

CH_color <- brewer.pal(4, "Spectral")[c(4)]
HH_color <- brewer.pal(4, "Spectral")[c(2)]
CD_color <- brewer.pal(4, "Spectral")[c(3)]
HD_color <- brewer.pal(4, "Spectral")[c(1)]
my_colors <- c(CH_color, HH_color, CD_color, HD_color)

CH_shp <- 15
HH_shp <- 19
CD_shp <- 22
HD_shp <- 21
CH_shp_box <- 22
HH_shp_box <- 21
my_shapes <- c(CH_shp, HH_shp, CD_shp, HD_shp)
my_shapes_box <- c(CH_shp_box, HH_shp_box, CD_shp, HD_shp)

my_labels <- c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")

Body Mass

Means Only F4

ggplot() + 
  #plot these first so they end up on the "bottom"
  geom_smooth(data = dat_no_rehab,
                  aes(x = day_n,
                y = mass_g, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means,
                  aes(x = day_n,
                y = mean_mass, 
                color = tmt,
                group = tmt,
                ymin = mean_mass-se_mass, 
                ymax = mean_mass+se_mass),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_mass, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(seq(9,12)),
                     labels = c(seq(9,12))) +
  xlab("Day") + 
  ylab("Body Mass (g)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 10.8, unit = "pt")
) -> mass_fig_min
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
mass_fig_min
## `geom_smooth()` using formula = 'y ~ x'

Ending Values F5

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = mass_g, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = mass_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = mass_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_y_continuous(limits = c(7,15),
                     breaks = c(seq(7,15, by = 2)),
                     labels = c(seq(7,15, by = 2))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  xlab("") + 
  
  annotate(geom = "text", x = 4, y = 14.4, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 15, label = "BC", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 14.6, label = "AC", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 14.8, label = "A", 
           size = 3) +
  
  ylab("Body Mass (g)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              3.4), "mm")
        ) -> mass_end_boxplot
mass_end_boxplot

mass_emmeans
##          tmt    emmean        SE       df  lower.CL upper.CL      response
## 1 Cool Humid 10.667647 0.2517552 46.72727 10.161103 11.17419 Body Mass (g)
## 2  Hot Humid  9.727273 0.2559428 47.91893  9.212643 10.24190 Body Mass (g)
## 3   Cool Dry 10.581818 0.2559914 47.91893 10.067090 11.09655 Body Mass (g)
## 4    Hot Dry  9.547059 0.2521465 46.16198  9.039562 10.05456 Body Mass (g)
mass_pairwise
##                 contrast    estimate        SE       df    t.ratio    p.value
## 1 Cool Humid - Hot Humid  0.94037433 0.3575715 126.2004  2.6298915 0.04674770
## 2  Cool Humid - Cool Dry  0.08582888 0.3585976 127.4556  0.2393459 0.99515381
## 3   Cool Humid - Hot Dry  1.12058824 0.3554223 126.8612  3.1528361 0.01074007
## 4   Hot Humid - Cool Dry -0.85454545 0.3615808 127.8215 -2.3633597 0.08951418
## 5    Hot Humid - Hot Dry  0.18021390 0.3583741 127.2163  0.5028653 0.95826222
## 6     Cool Dry - Hot Dry  1.03475936 0.3577985 126.4621  2.8920171 0.02307102
##        response
## 1 Body Mass (g)
## 2 Body Mass (g)
## 3 Body Mass (g)
## 4 Body Mass (g)
## 5 Body Mass (g)
## 6 Body Mass (g)

Hct

Ind + Means

ggplot() +
  geom_line(data = dat[complete.cases(dat$hematocrit_percent),],
            aes(x = day_n,
                y = hematocrit_percent, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means[complete.cases(means$mean_hct),],
            aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  geom_vline(xintercept = 9,
             linetype = "dashed",
             color = "darkgrey") +
  theme_classic() + 
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("") + 
  ylab("Hematocrit (%)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0.1, #top
                             0.1, #right
                             0.1, #bottom
                             0.41 #left
                             ), "cm")
        ) -> hct_fig
hct_fig
## Warning: Removed 12 rows containing missing values (`geom_point()`).

Means Only F5

ggplot() +
  geom_smooth(data = dat_no_rehab[complete.cases(dat_no_rehab$hematocrit_percent),],
                  aes(x = day_n,
                y = hematocrit_percent, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means[complete.cases(means$mean_hct),],
                  aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                group = tmt,
                ymin = mean_hct-se_hct, 
                ymax = mean_hct+se_hct),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means[complete.cases(means$mean_hct),],
            aes(x = day_n,
                y = mean_hct, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            alpha = 1, 
            fill = "white",
            size = 3) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(25, 30, 35, 40),
                     labels = c(25, 30, 35, 40),) +
  xlab("") + 
  ylab("Hematocrit (%)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 0, l = 10.8, unit = "pt")
        ) -> hct_fig_min
hct_fig_min
## `geom_smooth()` using formula = 'y ~ x'

Ending Values F5

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = hematocrit_percent, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = hct_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = hct_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  xlab("") + 
  scale_y_continuous(limits = c(10,55),
                     breaks = c(seq(10,50, by = 10)),
                     labels = c(seq(10,50, by = 10))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  
  annotate(geom = "text", x = 4, y = 47.6, label = "AB", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 44.6, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 54.6, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 54.6, label = "A", 
           size = 3) +
  
  ylab("Hematocrit (%)") + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              2.24), "mm")
        ) -> hct_end_boxplot
hct_end_boxplot
## Warning: Removed 3 rows containing missing values (`geom_point()`).

hct_emmeans
##          tmt   emmean       SE       df lower.CL upper.CL       response
## 1 Cool Humid 30.33352 1.808392 5.860156 25.88283 34.78422 Hematocrit (%)
## 2  Hot Humid 25.77326 1.830146 6.142261 21.32008 30.22645 Hematocrit (%)
## 3   Cool Dry 29.70335 1.815476 5.950456 25.25206 34.15464 Hematocrit (%)
## 4    Hot Dry 27.63092 1.815279 5.948653 23.17978 32.08205 Hematocrit (%)
hct_pairwise
##                 contrast   estimate       SE       df    t.ratio     p.value
## 1 Cool Humid - Hot Humid  4.5602630 1.272328 123.0398  3.5841894 0.002703439
## 2  Cool Humid - Cool Dry  0.6301769 1.254531 123.1123  0.5023206 0.958383983
## 3   Cool Humid - Hot Dry  2.7026074 1.253006 123.0665  2.1568990 0.141276348
## 4   Hot Humid - Cool Dry -3.9300861 1.288296 123.2298 -3.0506087 0.014665591
## 5    Hot Humid - Hot Dry -1.8576556 1.283706 123.1026 -1.4471033 0.472669199
## 6     Cool Dry - Hot Dry  2.0724305 1.262089 123.0607  1.6420637 0.359065321
##         response
## 1 Hematocrit (%)
## 2 Hematocrit (%)
## 3 Hematocrit (%)
## 4 Hematocrit (%)
## 5 Hematocrit (%)
## 6 Hematocrit (%)

Osml

Ind + Means

ggplot() + 
  geom_line(data = dat[complete.cases(dat$osmolality_mmol_kg_mean),],
            aes(x = day_n,
                y = osmolality_mmol_kg_mean, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means[complete.cases(means$mean_osml),],
            aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  geom_vline(xintercept = 9,
             linetype = "dashed",
             color = "darkgrey") +
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8, 10)) +
  ylim(300,450) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("") + 
  ylab("Plasma Osmolality (mmol/kg)") + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme_classic() +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0.6, #top
                             0.1, #right
                             0.1, #bottom
                             0.1 #left
                             ), "cm")
        ) -> osml_fig
osml_fig
## Warning: Removed 1 row containing missing values (`geom_line()`).
## Warning: Removed 12 rows containing missing values (`geom_point()`).

Means Only F4

ggplot() + 
  geom_smooth(data = dat_no_rehab[complete.cases(dat_no_rehab$osmolality_mmol_kg_mean),],
                  aes(x = day_n,
                y = osmolality_mmol_kg_mean, 
                color = tmt,
                group = tmt),
              method = "lm",
              se = F,
              size = 0.7) +
  geom_errorbar(data = means,
                  aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                group = tmt,
                ymin = mean_osml-se_osml, 
                ymax = mean_osml+se_osml),
                width = .1,
                alpha = 0.8) +
  geom_point(data = means,
            aes(x = day_n,
                y = mean_osml, 
                color = tmt,
                #fill = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(seq(340,380, by = 10)),
                     labels = c(seq(340,380, by = 10))) +
  xlab("") + 
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  #geom_vline(xintercept = 9,
   #          linetype = "dashed",
    #         color = "darkgrey") +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 0, l = 1, unit = "pt")
        ) -> osml_fig_min
osml_fig_min
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 12 rows containing missing values (`geom_point()`).

Stats! Check Pairwise Diffs ~ Time

Since Plasma osmolality has such a nonlinear trend, we need to test whether the elevated values in the middle of the experiment are significantly different than the values taken before and/or after.

# first make sub-dfs for each tmt group
HH <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Hu")
HD <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Hot Dr")
CH <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool H")
CD <- dat_no_rehab %>%
  dplyr::filter(substr(tmt, 1, 6) == "Cool D")

# next do pairwise tests for osml on the diff exp days, for each tmt group
pair_HH <- TukeyHSD(aov(data = HH, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_HD <- TukeyHSD(aov(data = HD, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_CH <- TukeyHSD(aov(data = CH, osmolality_mmol_kg_mean ~ day_factor)) #nonsig
pair_CD <- TukeyHSD(aov(data = CD, osmolality_mmol_kg_mean ~ day_factor)) #nonsig

# put into a df and export
osml_pairwise_df <- as.data.frame(pair_HD[[1]]) %>%
  rbind(as.data.frame(pair_HH[[1]])) %>%
  rbind(as.data.frame(pair_CD[[1]])) %>%
  rbind(as.data.frame(pair_CH[[1]])) %>%
  mutate(day_diff = paste("day", substr(rownames(.), 1, 3)),
         tmt = c(rep("Hot Dry",3),
                 rep("Hot Humid",3),
                 rep("Cool Dry",3),
                 rep("Cool Humid",3)),
         CI_95 = paste(round(lwr, digits = 2), round(upr, digits = 2), sep = ", "),
         diff = round(diff, digits = 2)) %>%
  dplyr::select(tmt, day_diff, diff, CI_95, p_adj = "p adj")

# save
write.csv(osml_pairwise_df, "./results_statistics/osmolality_pairwise_diffs.csv")

Nope, none of the differences between days within tmt groups are significantly different.

Ending Values F5

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = osmolality_mmol_kg_mean, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = osml_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = osml_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_y_continuous(limits = c(290,470),
                     breaks = c(seq(300,450, by = 50)),
                     labels = c(seq(300,450, by = 50))) +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  xlab("") + 
  
  annotate(geom = "text", x = 4, y = 427, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 452, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 437, label = "A", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 470, label = "A", 
           size = 3) +
  
  ylab(bquote('Osmolality (mmol '*kg^-1*')')) + 
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              0), "mm")
        ) -> osml_end_boxplot
osml_end_boxplot
## Warning: Removed 10 rows containing missing values (`geom_point()`).

osml_emmeans
##          tmt   emmean       SE       df lower.CL upper.CL
## 1 Cool Humid 350.3202 10.16423 4.874923 323.9893 376.6511
## 2  Hot Humid 354.0963 10.13690 4.824373 327.7506 380.4420
## 3   Cool Dry 349.1611 10.17830 4.903510 322.8413 375.4808
## 4    Hot Dry 362.0244 10.15727 4.863220 335.6918 388.3569
##                      response
## 1 Plasma Osmolality (mmol/kg)
## 2 Plasma Osmolality (mmol/kg)
## 3 Plasma Osmolality (mmol/kg)
## 4 Plasma Osmolality (mmol/kg)
osml_pairwise
##                 contrast   estimate       SE       df    t.ratio    p.value
## 1 Cool Humid - Hot Humid  -3.776107 5.031118 116.0115 -0.7505502 0.87626200
## 2  Cool Humid - Cool Dry   1.159115 5.128594 116.0589  0.2260104 0.99590771
## 3   Cool Humid - Hot Dry -11.704175 5.078569 116.0347 -2.3046207 0.10288296
## 4   Hot Humid - Cool Dry   4.935222 5.085161 116.0412  0.9705144 0.76645671
## 5    Hot Humid - Hot Dry  -7.928068 5.035558 116.0197 -1.5744172 0.39720576
## 6     Cool Dry - Hot Dry -12.863291 5.114607 116.0115 -2.5150105 0.06282759
##                      response
## 1 Plasma Osmolality (mmol/kg)
## 2 Plasma Osmolality (mmol/kg)
## 3 Plasma Osmolality (mmol/kg)
## 4 Plasma Osmolality (mmol/kg)
## 5 Plasma Osmolality (mmol/kg)
## 6 Plasma Osmolality (mmol/kg)

CEWL

Ind + Means

ggplot() + 
  geom_line(data = dat[complete.cases(dat$CEWL_g_m2h_mean),],
            aes(x = day_n,
                y = CEWL_g_m2h_mean, 
                color = tmt,
                group = individual_ID),
            alpha = 0.2) +
  geom_line(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                group = tmt),
            alpha = 1, 
            size = 1) +
  geom_point(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                shape = tmt),
            alpha = 1, 
            size = 5) +
  theme_classic() + 
  scale_shape_manual(values = c(15:18), name = "") +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_color_brewer(palette = "Set2", name = "") +
  xlab("Day") + 
  ylim(0,60) +
  ylab(bquote('CEWL (g/'*m^2*'h)')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 22),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 16),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 22),
        legend.text.align = 0,
        legend.position = "bottom"
        ) -> CEWL_fig
CEWL_fig

Means Only F2

ggplot() + 
  geom_errorbar(data = means[complete.cases(means$mean_CEWL),],
                  aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                group = tmt,
                ymin = mean_CEWL-se_CEWL, 
                ymax = mean_CEWL+se_CEWL),
                width = .1,
                #position=position_dodge(.1),
                alpha = 0.8) +
  geom_line(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                #linetype = tmt,
                group = tmt),
            alpha = 0.7, 
            size = .5) +
  geom_point(data = means[complete.cases(means$mean_CEWL),],
            aes(x = day_n,
                y = mean_CEWL, 
                color = tmt,
                shape = tmt),
            fill = "white",
            alpha = 1, 
            size = 3) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes, name = "",
                     labels = my_labels) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_color_manual(values = my_colors, name = "",
                     labels = my_labels) +
  scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
  scale_y_continuous(breaks = c(20, 25, 30, 35, 40),
                     limits = c(18,40)) +
  xlab("Day") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none"
        #legend.position = c(0.25,0.85)
        ) -> CEWL_fig_min
CEWL_fig_min

# use ggarrange so legend is centered
CEWL_fig_formatted <- ggarrange(CEWL_fig_min,
                                ncol = 1, nrow = 1,
                                common.legend = TRUE,
                                legend = "bottom")
# save
ggsave(filename = "experiment_CEWL_fig.pdf",
       plot = CEWL_fig_formatted,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 90)

delta ~ VPD F3

ggplot(data = dat_no_rehab_deltaCEWL) +
  geom_point(aes(x = VPD_mean_tmttrial,
                   y = delta_CEWL, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 2,
              alpha = 0.4
             ) +
  geom_smooth(aes(x = VPD_mean_tmttrial,
                   y = delta_CEWL), 
              se = F,
              formula = y ~ x,
              method = "lm",
              color = "black") +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_color_manual(values = my_colors, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_y_continuous(limits = c(-13,40),
                     breaks = c(seq(-10,40, by = 10)),
                     labels = c(seq(-10,40, by = 10))
                     ) +
  scale_x_continuous(limits = c(0,5),
                     breaks = c(0.6, 1.1, 2.5, 3.8),
                     labels = c("0.6\nCH", 
                                "1.1\nHH", 
                                "2.5\nCD", 
                                "3.8\nHD")) +
  xlab("Vapor Pressure Deficit (kPa)") + 
  ylab(expression(Delta ~ 'CEWL')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "bottom",
        legend.spacing.y = unit(0, "mm")
        #plot.margin = unit(c(0, #top
         #                     0, #right
          #                    0, #bottom
           #                   0), "mm")
        ) -> CEWL_VPD_fig
CEWL_VPD_fig

# code to export this figure alone
# NA bc now grouped with the next figure
# # use ggarrange so legend is centered
# CEWL_VPD_fig_formatted <- ggarrange(CEWL_VPD_fig,
#                                 ncol = 1, nrow = 1,
#                                 common.legend = TRUE,
#                                 legend = "bottom")
# # save
# ggsave(filename = "exp_CEWL_delta_VPD.pdf",
#        plot = CEWL_VPD_fig_formatted,
#        path = "./results_figures",
#        device = "pdf",
#        dpi = 600,
#        units = "mm",
#        width = 80, height = 90)

delta ~ WVP F3

ggplot(data = dat_no_rehab_deltaCEWL) +
  geom_point(aes(x = e_a_mean_tmttrial,
                   y = delta_CEWL, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 2,
              alpha = 0.4
             ) +
  geom_smooth(aes(x = e_a_mean_tmttrial,
                   y = delta_CEWL), 
              se = F,
              formula = y ~ x,
              method = "lm",
              color = "black") +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_fill_manual(values = my_colors, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_color_manual(values = my_colors, name = "",
                     labels = c("Cool Humid (CH)",
                               "Hot Humid (HH)",
                               "Cool Dry (CD)",
                               "Hot Dry (HD)")) +
  scale_y_continuous(limits = c(-13,40),
                     breaks = c(seq(-10,40, by = 10)),
                     labels = c(seq(-10,40, by = 10))
                     ) +
  scale_x_continuous(limits = c(0,5),
                     breaks = c(0.5, 2.0, 2.3, 4.5),
                     labels = c("0.5\nCD", 
                                "2.0  \nHD  ", 
                                "  2.3\n  CH", 
                                "4.5\nHH")) +
  xlab("Water Vapor Pressure (kPa)") + 
  ylab(expression(Delta ~ 'CEWL')) + 
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "bottom",
        legend.spacing.y = unit(0, "mm")
        #plot.margin = unit(c(0, #top
         #                     0, #right
          #                    0, #bottom
           #                   0), "mm")
        ) -> CEWL_WVP_fig
CEWL_WVP_fig

Ending Values F5

ggplot() +
  geom_jitter(data = end_vals,
              aes(x = tmt,
                   y = CEWL_g_m2h_mean, 
                   color = tmt,
                 fill = tmt,
                 shape = tmt), 
               size = 1,
              alpha = 0.4,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = CEWL_emmeans,
                aes(x = tmt,
                    y = emmean, 
                    color = tmt,
                    group = tmt,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 1) +
  geom_point(data = CEWL_emmeans, 
             aes(x = tmt,
                   y = emmean, 
                   #color = tmt,
                 shape = tmt,
                 fill = tmt), 
             color = "black",
               size = 4) +
  theme_classic() + 
  scale_shape_manual(values = my_shapes_box, name = "") +
  scale_fill_manual(values = my_colors, name = "") +
  scale_color_manual(values = my_colors, name = "") +
  scale_x_discrete(labels = c("Cool Humid\n0.6 kPa",
                               "Hot Humid\n1.1 kPa",
                               "Cool Dry\n2.5 kPa",
                               "Hot Dry\n3.8 kPa")) + 
  scale_y_continuous(limits = c(9,63),
                     breaks = c(seq(10,60, by = 10)),
                     labels = c(seq(10,60, by = 10))) +
  xlab("") + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) + 
  
  annotate(geom = "text", x = 4, y = 50, label = "C", 
           size = 3) +
  annotate(geom = "text", x = 2, y = 63, label = "B", 
           size = 3) +
  annotate(geom = "text", x = 3, y = 41, label = "C", 
           size = 3) +
  annotate(geom = "text", x = 1, y = 52, label = "A", 
           size = 3) +
  
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = unit(c(0, #top
                              0, #right
                              0, #bottom
                              0), "mm")
        ) -> CEWL_end_boxplot
CEWL_end_boxplot
## Warning: Removed 1 rows containing missing values (`geom_point()`).

CEWL_emmeans
##          tmt   emmean       SE       df lower.CL upper.CL     response
## 1 Cool Humid 30.36157 1.445712 10.68500 27.16810 33.55505 CEWL (g/m2h)
## 2  Hot Humid 36.76428 1.446643 10.69577 33.56916 39.95941 CEWL (g/m2h)
## 3   Cool Dry 23.72439 1.446548 10.69551 20.52946 26.91931 CEWL (g/m2h)
## 4    Hot Dry 24.61233 1.434910 10.37491 21.43074 27.79392 CEWL (g/m2h)
CEWL_pairwise
##                 contrast   estimate       SE       df    t.ratio      p.value
## 1 Cool Humid - Hot Humid -6.4027092 1.477747 125.0660 -4.3327496 1.740823e-04
## 2  Cool Humid - Cool Dry  6.6371867 1.482776 125.4155  4.4761900 9.849699e-05
## 3   Cool Humid - Hot Dry  5.7492451 1.468760 125.2015  3.9143543 8.436015e-04
## 4   Hot Humid - Cool Dry 13.0398959 1.482912 125.4316  8.7934385 9.026113e-14
## 5    Hot Humid - Hot Dry 12.1519543 1.469693 125.2717  8.2683648 1.057265e-12
## 6     Cool Dry - Hot Dry -0.8879416 1.467098 125.0811 -0.6052366 9.302438e-01
##       response
## 1 CEWL (g/m2h)
## 2 CEWL (g/m2h)
## 3 CEWL (g/m2h)
## 4 CEWL (g/m2h)
## 5 CEWL (g/m2h)
## 6 CEWL (g/m2h)

Exp CEWL ~ Osml

end_vals_CEWL_osml <- dat %>%
  dplyr::filter(day_n == 8) %>%
  dplyr::filter(complete.cases(CEWL_g_m2h_mean, osmolality_mmol_kg_mean))

ggplot(end_vals_CEWL_osml) + 
  aes(x = osmolality_mmol_kg_mean,
      y = CEWL_g_m2h_mean,
      color = tmt,
      shape = tmt) +
  geom_point(size = 1.5, 
             alpha = 0.8) + 
  stat_smooth(formula = y ~ x, 
              method = "lm", 
              se = F, 
              size = 1, 
              alpha = 0.9) + 
  theme_classic() + 
  xlab(bquote('Osmolality (mmol '*kg^-1*')')) + 
  ylab(bquote('CEWL (g '*m^-2*' '*h^-1*')')) + 
  scale_shape_manual(values = c(21,19, 22,15), name = "") +
  scale_fill_brewer(palette = "Spectral", name = "") +
  scale_color_brewer(palette = "Spectral", name = "") +
  scale_x_continuous(breaks = c(300, 350, 400, 450)) +
  scale_y_continuous(breaks = c(20, 30, 40, 50),
                     limits = c(12,57)) +
  guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.position = "bottom"
        ) -> exp_end_CEWL_osml_fig
exp_end_CEWL_osml_fig

ggsave(filename = "exp_CEWL_osml_fig.pdf",
       plot = exp_end_CEWL_osml_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 90)

Multi-Figures

hydration ~ time F4

ggarrange(osml_fig_min, 
          hct_fig_min, 
          mass_fig_min,
          ncol = 1, nrow = 3,
          labels = c("A", "B", "C"),
          common.legend = TRUE,
          legend = "bottom"
          ) -> experiment_multi_fig
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 12 rows containing missing values (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 12 rows containing missing values (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#experiment_multi_fig

# export figure
ggsave(filename = "experiment_multi_fig.pdf",
       plot = experiment_multi_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 210)

end values F5

ggarrange(CEWL_end_boxplot, 
          osml_end_boxplot, 
          hct_end_boxplot,
          mass_end_boxplot, 
          ncol = 2, nrow = 2,
          labels = c("A", "B", "C", "D"),
          widths = c(2, 2.045), heights = c(2, 2),
          common.legend = FALSE
          ) -> ending_values_multi_fig
## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Warning: Removed 10 rows containing missing values (`geom_point()`).
## Warning: Removed 3 rows containing missing values (`geom_point()`).
#ending_values_multi_fig

ggsave(filename = "exp_end_val_multi_fig.pdf",
       plot = ending_values_multi_fig,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 180, height = 150)

VPD/WVP F3

# use ggarrange so legend is centered
CEWL_VPD_fig_formatted <- ggarrange(CEWL_VPD_fig,
                                    CEWL_WVP_fig,
                                ncol = 1, nrow = 2,
                                common.legend = TRUE,
                                legend = "bottom",
                                labels = c("A", "B"))
# save
ggsave(filename = "exp_CEWL_delta_VPD_WVP.pdf",
       plot = CEWL_VPD_fig_formatted,
       path = "./results_figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 150)